Locosselli et al. PNAS modeling analysis

# libraries
library(MCMCglmm); library(coda)
library(parallel)
library(phytools); library(brranching); library(ape)
library(dplyr)

Dataset

df <- read.csv("data/dataset.csv")
str(df)
## 'data.frame':    763 obs. of  18 variables:
##  $ zone        : Factor w/ 3 levels "Sub-tropics",..: 3 3 3 3 3 3 3 3 3 3 ...
##  $ ID          : int  1 2 3 4 5 6 7 8 9 10 ...
##  $ fam         : Factor w/ 53 levels "Anacardiaceae",..: 36 30 25 25 25 25 21 27 30 21 ...
##  $ taxagroup   : Factor w/ 2 levels "Conifer","Eudicot": 1 2 2 2 2 2 2 2 2 2 ...
##  $ species     : Factor w/ 281 levels "Abies delavayi",..: 5 267 260 260 260 260 26 44 71 73 ...
##  $ longevity   : int  311 186 221 523 262 219 243 427 308 123 ...
##  $ growth      : num  1.94 5.52 4 2.62 NA ...
##  $ densityCHAVE: num  0.43 0.56 0.601 0.601 0.601 0.601 0.524 0.638 0.447 0.504 ...
##  $ human.influ : int  35 39 56 14 18 22 12 15 12 15 ...
##  $ tseas       : int  1340 1669 1941 1043 1325 1036 1095 685 1095 685 ...
##  $ moistwetq   : num  1.18 1.38 2.47 1.96 2.19 1.59 1.32 1.36 1.32 1.36 ...
##  $ low.w.rad   : num  145 145 136 144 143 ...
##  $ temp        : num  24.5 22.1 27 22.6 27.7 21.6 25.1 26.5 25.1 26.5 ...
##  $ moistdq     : num  0.18 0.19 0.09 0.4 0.25 0.52 0.44 0.34 0.44 0.34 ...
##  $ cloud       : int  302 307 596 207 329 179 157 210 157 210 ...
##  $ soil_class  : Factor w/ 26 levels "Acrisols","Alisols",..: 13 16 9 19 16 16 9 8 9 8 ...
##  $ alt.group   : Factor w/ 2 levels "High altitude",..: 1 2 2 2 2 2 2 2 2 2 ...
##  $ phylo       : Factor w/ 281 levels "Abies_delavayi",..: 5 267 260 260 260 260 26 44 71 73 ...

Columns information/metadata:

  • zone: global climatic zone
  • ID: identity of the observation data
  • fam: Family name
  • taxagroup: taxonomical group, Conifer or Eudicot
  • species: Species name
  • longevity: recorde maximum longevity
  • growth: recorder growth
  • densityCHAVE: wood density
  • human.influ: Human influence index
  • tseas: temperature seasonality
  • moistwetq: water-balance soil moisture index of the wettest quarter of the year
  • low.w.rad: lowest weekly radiation
  • temp: mean annual temperature
  • moistdq:moisture index of the driest quarter
  • cloud: cloud cover
  • soil_class: soil classification
  • alt.group: separation of the populations in “High altitude” and “Lowland”
  • phylo: Species names used to match the phylogenetic tree

Standardizing numeric variables:

data <- df %>% mutate_at(vars(8:15), scale)

Phylogeny

The phylogenetic tree was constructed using the Phylomatic software (Webb et al. 2008) based on a revised vascular megatree (Gastauer & Meira-Neto 2016). We estimated the ancestral values of longevity and growth at internal nodes using maximum likelihood under the assumption of Brownian motion for trait evolution (Revel 2013). Auxiliary packages used: phytools, brranching and ape.

filo <- read.tree("data/phylogenetic_tree.new")

Inverse phylogenetic matrix to be used in the model

inv.phylo <- inverseA(filo, nodes = "TIPS", scale = FALSE)

Model 1

Code of the bivariate phylogenetic linear mixed model with the packages MCMCglmm and coda.

The model takes long yours (or days) to run. The final version was stored to be loaded:

load("model_output/modelResults.Rdata")

Prior specification for the random effects and residuals

prior <- list(R = list(V = diag(2), nu = 1.002),
               G = list(G1 = list(V = diag(2), nu = 1.002), # phylogeny
                        G2 = list(V = diag(2), nu = 1.002), # species
                        G3 = list(V = diag(2), nu = 1.002))) # soil_class

Model

mod <- MCMCglmm(cbind(log(longevity), log(growth)) ~ trait
        + trait:human.influ
        + trait:densityCHAVE 
        + trait:moistwetq + trait:temp + trait:moistdq
        + trait:tseas + trait:low.w.rad + trait:cloud - 1,
        random = ~ us(trait):phylo + us(trait):species + idh(trait):soil_class,
        rcov = ~ us(trait):units,
        family = c("gaussian", "gaussian"), 
        ginverse = list(phylo = inv.phylo$Ainv),
        prior = prior, data = data, pr = TRUE,
        nitt = 1000000, burnin = 10000, thin = 200)

Model summary

summary(mod)
## 
##  Iterations = 10001:999801
##  Thinning interval  = 200
##  Sample size  = 4950 
## 
##  DIC: 1512.182 
## 
##  G-structure:  ~us(trait):phylo
## 
##                                     post.mean l-95% CI u-95% CI eff.samp
## traitlongevity:traitlongevity.phylo    0.4993  0.19931   0.7953     4950
## traitgrowth:traitlongevity.phylo      -0.1955 -0.39760  -0.0275     4950
## traitlongevity:traitgrowth.phylo      -0.1955 -0.39760  -0.0275     4950
## traitgrowth:traitgrowth.phylo          0.2062  0.06251   0.3730     4950
## 
##                ~us(trait):species
## 
##                                       post.mean l-95% CI u-95% CI eff.samp
## traitlongevity:traitlongevity.species   0.10901  0.04823  0.17329     4950
## traitgrowth:traitlongevity.species     -0.02501 -0.06778  0.01552     4950
## traitlongevity:traitgrowth.species     -0.02501 -0.06778  0.01552     4950
## traitgrowth:traitgrowth.species         0.09071  0.04905  0.13286     4739
## 
##                ~idh(trait):soil_class
## 
##                           post.mean l-95% CI u-95% CI eff.samp
## traitlongevity.soil_class   0.06294  0.02436  0.10914     4950
## traitgrowth.soil_class      0.05488  0.02245  0.09276     5420
## 
##  R-structure:  ~us(trait):units
## 
##                                     post.mean l-95% CI u-95% CI eff.samp
## traitlongevity:traitlongevity.units   0.18907  0.16606  0.21423     4950
## traitgrowth:traitlongevity.units     -0.05821 -0.07277 -0.04211     4950
## traitlongevity:traitgrowth.units     -0.05821 -0.07277 -0.04211     4950
## traitgrowth:traitgrowth.units         0.11571  0.09917  0.13324     4950
## 
##  Location effects: cbind(log(longevity), log(growth)) ~ trait + trait:human.influ + trait:densityCHAVE + trait:moistwetq + trait:temp + trait:moistdq + trait:tseas + trait:low.w.rad + trait:cloud - 1 
## 
##                             post.mean  l-95% CI  u-95% CI eff.samp    pMCMC    
## traitlongevity               5.172109  4.499176  5.821662     4950  < 2e-04 ***
## traitgrowth                  1.129257  0.703709  1.586139     5153  < 2e-04 ***
## traitlongevity:human.influ  -0.105172 -0.149308 -0.061976     4950  < 2e-04 ***
## traitgrowth:human.influ     -0.023794 -0.064216  0.013438     4310 0.220606    
## traitlongevity:densityCHAVE -0.046698 -0.121774  0.031116     4950 0.242020    
## traitgrowth:densityCHAVE    -0.121748 -0.185407 -0.057769     4950  < 2e-04 ***
## traitlongevity:moistwetq     0.037342 -0.024483  0.100625     4950 0.258182    
## traitgrowth:moistwetq        0.038007 -0.015371  0.091839     4950 0.174141    
## traitlongevity:temp          0.031938 -0.052077  0.114523     4950 0.454545    
## traitgrowth:temp             0.107223  0.033680  0.187452     5202 0.005253 ** 
## traitlongevity:moistdq       0.148818  0.070339  0.226906     4714 0.000404 ***
## traitgrowth:moistdq         -0.035666 -0.104497  0.028710     5462 0.294545    
## traitlongevity:tseas         0.049887 -0.037768  0.135077     4950 0.267879    
## traitgrowth:tseas           -0.073050 -0.149220  0.008314     4950 0.074747 .  
## traitlongevity:low.w.rad    -0.083994 -0.148693 -0.021337     4950 0.012121 *  
## traitgrowth:low.w.rad        0.041644 -0.011022  0.094039     4950 0.113939    
## traitlongevity:cloud        -0.013544 -0.075064  0.054227     4950 0.685657    
## traitgrowth:cloud            0.034173 -0.026544  0.088690     5647 0.239192    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Gelman-Rubin criteria of convergence (4 chains)

4 chains of the model previously ran to make the diagnose of the convergence of the chain:

gelman.diag(mcov4[,1:18])
## Potential scale reduction factors:
## 
##                             Point est. Upper C.I.
## traitlongevity                       1          1
## traitgrowth                          1          1
## traitlongevity:human.influ           1          1
## traitgrowth:human.influ              1          1
## traitlongevity:densityCHAVE          1          1
## traitgrowth:densityCHAVE             1          1
## traitlongevity:moistwetq             1          1
## traitgrowth:moistwetq                1          1
## traitlongevity:temp                  1          1
## traitgrowth:temp                     1          1
## traitlongevity:moistdq               1          1
## traitgrowth:moistdq                  1          1
## traitlongevity:tseas                 1          1
## traitgrowth:tseas                    1          1
## traitlongevity:low.w.rad             1          1
## traitgrowth:low.w.rad                1          1
## traitlongevity:cloud                 1          1
## traitgrowth:cloud                    1          1
## 
## Multivariate psrf
## 
## 1

Trace and density plots for fixed effects

par(mfrow=c(4,2), mar=c(2, 1, 1, 1))
plot(mcov4[,1:18], ask=F, auto.layout=F)

Trace and density plots for random effects

par(mfrow=c(4,2), mar=c(2, 1, 1, 1))
plot(mcov4.var, ask=F, auto.layout=F)

Model 2 lowland

Subsetting data:

datalow <- data[data$alt.group == "Lowland" & df$temp >=20, ]

Prior specification for the random effects and residuals

prior <- list(R = list(V = diag(2), nu = 1.002),
               G = list(G1 = list(V = diag(2), nu = 1.002), # phylogeny
                        G2 = list(V = diag(2), nu = 1.002), # species
                        G3 = list(V = diag(2), nu = 1.002))) # soil_class

Model

mod.low20 <- MCMCglmm(cbind(log(longevity), log(growth)) ~ trait +
                      + trait:temp + trait:temp2 + trait:temp3 - 1,
                      random = ~ us(trait):phylo + us(trait):species + idh(trait):soil_class,
                      rcov = ~ us(trait):units,
                      family = c("gaussian", "gaussian"), ginverse = list(phylo = inv.phylo$Ainv),
                      prior = prior, data = datalow, pr = TRUE,
                      nitt = 100000, burnin = 10000, thin = 200)

Model’s results previoulsy ran:

load("model_output/modelResult_lowland.Rdata")

Model summary

summary(mod.low20)
## 
##  Iterations = 10001:99801
##  Thinning interval  = 200
##  Sample size  = 450 
## 
##  DIC: 1015.704 
## 
##  G-structure:  ~us(trait):phylo
## 
##                                     post.mean l-95% CI u-95% CI eff.samp
## traitlongevity:traitlongevity.phylo    0.8533  0.32932  1.39428      450
## traitgrowth:traitlongevity.phylo      -0.2310 -0.51216  0.01856      450
## traitlongevity:traitgrowth.phylo      -0.2310 -0.51216  0.01856      450
## traitgrowth:traitgrowth.phylo          0.2103  0.04624  0.39553      450
## 
##                ~us(trait):species
## 
##                                       post.mean l-95% CI u-95% CI eff.samp
## traitlongevity:traitlongevity.species    0.1316  0.04960  0.23736      450
## traitgrowth:traitlongevity.species      -0.0235 -0.08112  0.03525      450
## traitlongevity:traitgrowth.species      -0.0235 -0.08112  0.03525      450
## traitgrowth:traitgrowth.species          0.1043  0.05968  0.16160      450
## 
##                ~idh(trait):soil_class
## 
##                           post.mean l-95% CI u-95% CI eff.samp
## traitlongevity.soil_class   0.07081  0.02655  0.14000      450
## traitgrowth.soil_class      0.05778  0.02368  0.09976      450
## 
##  R-structure:  ~us(trait):units
## 
##                                     post.mean l-95% CI u-95% CI eff.samp
## traitlongevity:traitlongevity.units   0.20939  0.17288   0.2472    450.0
## traitgrowth:traitlongevity.units     -0.06491 -0.09168  -0.0421    486.5
## traitlongevity:traitgrowth.units     -0.06491 -0.09168  -0.0421    486.5
## traitgrowth:traitgrowth.units         0.12699  0.10400   0.1507    679.8
## 
##  Location effects: cbind(log(longevity), log(growth)) ~ trait + +trait:temp + trait:temp2 + trait:temp3 - 1 
## 
##                      post.mean l-95% CI u-95% CI eff.samp  pMCMC   
## traitlongevity          5.2951   4.2999   6.1978      450 <0.002 **
## traitgrowth             1.2774   0.6980   1.7618      450 <0.002 **
## traitlongevity:temp    -2.4071  -4.5952  -0.2509      450 0.0267 * 
## traitgrowth:temp       -0.4205  -1.7737   1.4981      450 0.5867   
## traitlongevity:temp2    5.3139   1.8497   8.8137      450 <0.002 **
## traitgrowth:temp2       1.1204  -1.7590   3.8290      450 0.4400   
## traitlongevity:temp3   -3.1367  -4.7923  -1.4586      450 <0.002 **
## traitgrowth:temp3      -0.5973  -1.9841   0.7990      450 0.3956   
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Trace and density plots for fixed effects

par(mfrow=c(4,2), mar=c(2, 1, 1, 1))
plot(mod.low20$Sol[,1:8], ask=F, auto.layout=F)

Trace and density plots for random effects

par(mfrow=c(4,2), mar=c(2, 1, 1, 1))
plot(mod.low20$VCV, ask=F, auto.layout=F)

Session Information

sessionInfo()
## R version 3.6.2 (2019-12-12)
## Platform: x86_64-apple-darwin15.6.0 (64-bit)
## Running under: macOS Catalina 10.15.7
## 
## Matrix products: default
## BLAS:   /Library/Frameworks/R.framework/Versions/3.6/Resources/lib/libRblas.0.dylib
## LAPACK: /Library/Frameworks/R.framework/Versions/3.6/Resources/lib/libRlapack.dylib
## 
## locale:
## [1] pt_BR.UTF-8/pt_BR.UTF-8/pt_BR.UTF-8/C/pt_BR.UTF-8/pt_BR.UTF-8
## 
## attached base packages:
## [1] parallel  stats     graphics  grDevices utils     datasets  methods   base     
## 
## other attached packages:
## [1] dplyr_1.0.1      brranching_0.5.0 phytools_0.6-99  maps_3.3.0       MCMCglmm_2.29   
## [6] ape_5.3          coda_0.19-3      Matrix_1.2-18    knitr_1.23      
## 
## loaded via a namespace (and not attached):
##  [1] jsonlite_1.6.1          foreach_1.4.4           bold_0.9.0             
##  [4] gtools_3.8.1            shiny_1.3.2             assertthat_0.2.1       
##  [7] expm_0.999-4            highr_0.8               animation_2.6          
## [10] tensorA_0.36.1          taxize_0.9.8            yaml_2.2.0             
## [13] pillar_1.4.6            numDeriv_2016.8-1.1     lattice_0.20-38        
## [16] glue_1.4.2              quadprog_1.5-7          phangorn_2.5.5         
## [19] uuid_0.1-2              digest_0.6.25           promises_1.0.1         
## [22] htmltools_0.3.6         httpuv_1.5.1            plyr_1.8.4             
## [25] pkgconfig_2.0.3         httpcode_0.2.0          questionr_0.7.0        
## [28] bookdown_0.11           purrr_0.3.4             xtable_1.8-4           
## [31] corpcor_1.6.9           later_0.8.0             phylocomr_0.1.4        
## [34] cubature_2.0.3          tibble_3.0.3            combinat_0.0-8         
## [37] generics_0.0.2          ellipsis_0.3.1          cli_2.0.2              
## [40] mnormt_1.5-5            magrittr_1.5            crayon_1.3.4           
## [43] mime_0.7                evaluate_0.14           fansi_0.4.1            
## [46] nlme_3.1-142            MASS_7.3-51.4           xml2_1.3.2             
## [49] tools_3.6.2             data.table_1.12.8       lifecycle_0.2.0        
## [52] stringr_1.4.0           plotrix_3.7-8           compiler_3.6.2         
## [55] rlang_0.4.7             clusterGeneration_1.3.4 grid_3.6.2             
## [58] conditionz_0.1.0        iterators_1.0.10        rstudioapi_0.11        
## [61] igraph_1.2.4.2          miniUI_0.1.1.1          rmarkdown_1.13         
## [64] codetools_0.2-16        reshape_0.8.8           curl_3.3               
## [67] reshape2_1.4.3          R6_2.4.1                zoo_1.8-6              
## [70] fastmatch_1.1-0         stringi_1.4.3           rmdformats_0.3.3       
## [73] crul_0.8.4              Rcpp_1.0.4              vctrs_0.3.4            
## [76] tidyselect_1.1.0        scatterplot3d_0.3-41    xfun_0.7

References

Webb C.O., D. D. Akerly, S. W. Kembel, Phylocom: software for the analysis of phylogenetic community structure and character evolution. Bioinformatics, 24, 2098-2100. http://dx.doi.org/10.1093/bioinformatics/btn358. PMid:18678590 (2008).

Gastauer,M., J. A. A. Meira-Neto, An enhanced calibration of a recently released megatree for the analysis of phylogenetic diversity. Brazilian Journal of Biology, 76(3), 619-628 (2016).

Revell, L.J. Two new graphical methods for mapping trait evolution on phylogenies. Methods in Ecology and Evolution 4, 754–759 (2013).

Melina Leite

17 de Outubro de 2020